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Abstract. A rigorous convergence analysis of the Strang splitting algorithm with a discontinuous 
Galerkin approximation in space for the Vlasov-Poisson equations is provided. It is shown that under 
suitable assumptions the error is of order O (r^ + /i'' + h"^ /t^, where r is the size of a time step, h is 
the cell size, and q the order of the discontinuous Galerkin approximation. In order to investigate the 
recurrence phenomena for approximations of higher order as well as to compare the algorithm with 
numerical results already available in the literature a number of numerical simulations are performed. 
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1. Introduction. In astro- and plasma physics the behavior of a colhsionless 
plasma is modeled by the Vlasov equation (see e.g. [2]) 

dtf{t, x,v)+v V/(t, x,v) + F- V^f{t, X, v) = 0, (1.1) 

a kinetic model that in certain applications is also called the colhsionless Boltzmann 
equation. It is posed in a 3 + 3 dimensional phase space, where x denotes the position 
and V the velocity. The density function / is the sought-after particle-probability 
distribution, and the (force) term F describes the interaction of the plasma with the 
electromagnetic field. 

In this paper we will study the convergence properties of a full discretization of 
the so called Vlasov-Poisson equations, where the force term 

F = -V^ 

is the gradient of the self-consistent electric potential 0. This simplified model is used 
in various applications, e.g. in the context of Landau damping. 



For the numerical solution of (1.1), various methods have been considered in the 
literature, for example particle methods and in particular the particle-in-cell method, 
see [lOl [n] [T^ . Another prevalent approach consists in employing splitting methods, 
first proposed in the context of the Vlasov-Poisson equations by [7] and later extended 
to the full Vlasov-Maxwell equations in [T7]. Both papers use second-order Strang 
splitting. In the seminal paper [TS], the convergence properties of Strang-splitting for 
evolution equations were analyzed with the help of the variation-of-constants formula. 
This approach was recently extended to Vlasov- type equations in [Qj. In jSHUIjQ] semi- 
Lagrangian methods are combined with Strang splitting. Convergence is shown in the 
case of the 1-1-1 dimensional Vlasov-Poisson equations, and in [S| for a special case 
of the one-dimensional Vlasov-Maxwell equation. In these papers usually Hermite or 
spline interpolation is employed. 

On the other hand, discontinuous Galerkin approximations in space have been 
studied for the Vlasov-Poisson equations as well. In [13] and [14] the weak version of 
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the Vlasov-Poisson equations is discretized by a discontinuous Galerkin scheme and 
integrated in time by Runge-Kutta methods. In |^ a higher order semi-Lagrangian 
method in time is combined with a discontinuous Galerkin approximation in space. 
However, no convergence analysis is given. A direct Strang splitting scheme with 
a discontinuous Galerkin approximation is implemented in [17]. Since only a single 
value per cell is advanced in time this leads to a Van Leer scheme. In the before 
mentioned paper a numerical study of this scheme is conducted. 

In this paper, we will extend the analysis done in [5] for the Vlasov-Poisson 
equations in 1+1 dimensions to the fully discretized case. More precisely, we will 
show that the (direct) Strang splitting scheme combined with a discontinuous Galerkin 
discretization in space is convergent. Our main result is stated in Theorem |2.9| below . 
In addition, we will discuss some numerical aspects of our discretization in section [3j 

2. Vlasov-Poisson equations in 1+1 dimensions. In this section we per- 
form the convergence analysis of Strang splitting in time with a discontinuous Galerkin 
approximation in space for the Vlasov-Poisson equations in 1+1 dimensions. To that 
end we first describe the setting as well as give the necessary regularity results (sections 
2.1|to 2.3). We then describe the time (section 2.4) and space discretization (sections 



2.5 



and 2.6). In section 2.7 we will extend a commonly employed approximation re- 



sult from C^"*"^ functions to piecewise polynomials with a small jump discontinuity. 
This extension is crucial to show consistency (which is done in section 2.8). Finally, 



convergence is established in section 2.9 



2.1. Setting. We will consider the Vlasov-Poisson equations in 1+1 dimensions. 



I.e. 



dtf{t, X, v) = -vdxf{t, X, v) - E{f{t, •, •), x)dyf{t, x, v) 



d^E{f{t,;-),x) 



f(t,x,v) dv-l 



(2.1) 



f{0,x,v) = foix,v) 



with periodic boundary conditions in space. The domain of interest is given by 
{t,x,v) S [0,T] X [0, L] X M. The periodic boundary conditions imply 

Va; e M: f{t,x,v) = f[t,x + L,v). 

It is physically reasonable to assume at least an algebraic decay of /o in the velocity 
direction. Thus, we can approximate (to arbitrary precision) /o by an initial value 
with compact support. As will be apparent in the next section it is unnecessary to 
impose boundary conditions in the velocity direction for initial values with compact 
support. This is due to the fact that for such an initial value the solution will continue 



to have compact support for all finite time intervals [0,T] (see Theorem 2.1). 

For most of this presentation it will be more convenient to work with the following 
abstract initial value problem 



dtf{t) = {A + B)f{t) 
/(O) = /o, 



(2.2) 



where we assume that A is an (unbounded) linear operator. In addition, we assume 
that B can be written in the form Bf = B{f)f, where B{f) is an (unbounded) linear 
operator. For the Vlasov-Poisson equation in 1+1 dimensions the obvious choice is 
Af = -vd^f and Bf ^ ~E{f{t, •, •), x)dj. 
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In 1 + 1 dimensions an explicit representation of the electric field is given by the 
following formula 



E{f{t,;-),x)= I K{x,y)\ I f{t,y,v)dv~l]dy, 



K{x,y) 




(2.3) 



where E is chosen to have zero integral mean (electrostatic condition) . This represen- 
tation allows us to get a simple estimate of the electric field in terms of the probability 
density /. Note, however, that all the estimates employed in this paper could just as 
well be derived from potential theory. In fact, this approach is preferred as soon as 
one considers more than a single dimension in the space direction. 

2.2. Definitions and notation. The purpose of this section is to introduce the 
notations and mathematical spaces necessary for giving existence, uniqueness, and 
regularity results as well as to conduct the estimates necessary for showing consistency 
and convergence. 

First we introduce some notations that will be employed in the subsequent anal- 
ysis. Suppose that the differential equation g' — G{g) has (for a given initial value) 
a unique solution. In this case we denote the solution at time t with initial value 
5(^0) = 50 by g{t) = Ea{t — to, .9o)- In addition, we will use || • || to denote the infinity 
norm and || • ||p to denote the norm on [0,L] x M. 

For estimating the errors in space and velocity we will use the Banach space 
L°°([0, L] X [— Umax, Winax])- Notc that consistcucy bounds in the physically more 
reasonable norm are a direct consequence of the bounds we derive in the infinity 
norm. The situation is more involved in the case of stability (this is discussed in 
section 



2.91 



For our convergence analysis we need some regularity of the solution. To that 
end, we introduce the following spaces of continuously differentiable functions 

C^r,c {9 ^ C™(M^M), yx,v: {g{x + L,v) = g{x,v)) A (supp5i(a;, •) compact)} , 
C™, {g e C"\R,R), Vx: g{x + L) ^ g{x)} . 

Equipped with the norm of uniform convergence of all derivatives up to order to, C™j. ^ 
and Cpcr are Banach spaces. 

We also have to consider spaces that involve time. To that end let us define for 
any subspace Z C C""(M'',M) the space 

C"(0, T; Z) 1 5 e C™([0, T], C°), {g{t) e Z) A ( sup \\g{t)\\^ < 00) I . 

i *e[0:T] J 

Below, we will either take the choice Z = C™j. ^ or Z = C™^. It should be noted that 
functions in C"^(0,T; Z) possess spatial derivatives up to order to that are uniformly 
bounded in t e [0,T]. 

2.3. Existence, uniqueness, and regularity. In this section we recall the 
existence, uniqueness, and regularity results of the Vlasov-Poisson equations in 1+1 
dimensions. The theorem is stated with a slightly different notation in ^ and [3^. 

Theorem 2.1. Assume that /o £ C^^.c non-negative, then f e ^"(0, T;C™r_(.) 
and E{f{t,-,-),x) as a function of{t,x) lies m C™(0,T;C™,.). In addition, we can 
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find a number Q(T) such that for all t € [0, T] and a; € M it holds that supp /(t, a;, •) C 

[-g(r),Q(T)]. 

Proof. A proof can be found in [T2] • □ 

We also need a regularity result for the electric field that does not directly result 
from a solution of the Vlasov-Poisson equations, but from some generic function / 
(e.g., an / computed from an application of a splitting operator to /o). 

Corollary 2.2. For f e C;;^,.^ it holds that E{f, ■) e C"^ 

Proof. The result follows from the proof of Theorem |2.1| In addition, in the 
1+1 dimensional case it can also be followed from the exact representation of the 



electromagnetic field that is given in equation (2.3 1. □ 

It should also be noted that due to the proof of Theorem |2.1[ the regularity results 
given can be extended to the differential equ ations generated by B and B{g) (for any 
sufficiently regular g). Thus, Theorem 2.1 remains valid if Esit^fo) or e*'^'^-'/o is 



substituted for f{t). 

2.4. Time discretization. We use Strang-splitting for the time discretization 



of (2.2 1. This results in the scheme 

fk+i = Skfk, (2.4a) 

where fk is the numerical approximation to f{t) at time t = kr with step size r. The 
splitting operator Sk is the composition of three abstract operators 

5fc = (2.4b) 

where 

= e5^, ^ = e-^W-^+v^) (2.4c) 

with fk+i/2 — e^^'^^'^^ei'^fk. The choice of fk+1/2 is such as to retain second order 
in the non-linear case while still only advection problems have to be solved in the 
numerical approximation (for more details see e.g. [9]). Note that since e^^'^^^'^ can 
be represented by a translation in the velocity direction only (which has no effect on 
the computation of the electric field) we can use here 

/fe+1/2 - S^^^fu. (2.4d) 

This is convenient as the computation of S^^^ fk incurs no performance overhead in 
the actual computation. 

2.5. Space discretization. We proceed in two steps. First, we introduce a 
cutoff in the velocity direction, i.e. we fix Wmax and consider the problem on the 
domain [0,L] x [— fmax, Wmax]- Note that for an initial value with compact support 
with respect to velocity and a sufficiently large Vmax this is still exact. 

Second, we introduce a discontinuous Galerkin approximation in both the space 
and velocity direction. For simplicity, we consider a uniform rectangular grid. In this 
case, the cell boundaries are given by the following relations 

Xi = ihx, <i < Nx, 

Vj = jh^ - Umax, < .j < Ny. 

Within each ceU, i.e. a square Rij = [ih^, [i + l)/ix] x [jK - Wmax, (j + l)/iti - Wmax], 
< i < Nx, < j < Ny, we perform an orthogonal projection with respect to the 
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basis of Legendre polynomials of degree at most i in x and v. To be more precise, 
suppose that g e L^([0,i] x [— I'max, t'max]); then the operator P is defined such 
that Pg restricted to Rij for all i, j is the (unique) polynomial that results from 
the projection of g onto the {£ + + 1) dimensional subspace generated by the 
(appropriately translated and scaled) Legendre polynomials up to degree £. It is well 
known that this projection operator is given by 

PgK = E E bLH'\^)P^Hv) (2.5a) 

k=0 m=0 

with coefficients 

^rn- ^'^'^f"^'^ I P!^\x)Pj^'\v)g{a:,v)dia:,v). (2.5b) 



The translated and scaled Legendre polynomials are here defined as 

where pi denote the Legendre polynomials with the standard normalization, i.e. 

2 



-1 



pi{x)pj{x)dx = i^j^Sij. 



It should be emphasized that the projection in a single cell is independent from 
the projection in any other cell. As this is not true for Hcrmitc or spline interpolation 
it gives the discontinuous Galerkin scheme a computational advantage (see |17j and 
section 



2.6). 



Now we have to introduce an approximation to the abstract splitting operator 



(2.4b) that takes the space discretization into account. We use the decomposition 

Sk^S^^^S^^^S^^K (2.6a) 

where 

SiA) ^ pg{A)^ g{B) ^ p^rB(h^,/,) (2.6b) 

with 

/fe+i/2 = S^^^h- (2.6c) 
The fully discrete scheme then reads 

h+i = Skfk, /o - Pfo- (2.6d) 
Note that /j. represents the full approximation in time and space at time t^. 

2.6. Translation and projection. The principle algorithm has already been 
laid out in sections 2.4 and 2.5 However, the description given so far is certainly not 
sufficient as the straightforward implementation (first computing an exact solution 
and then projection onto a finite dimensional subspace) is clearly not a viable option. 
Thus, the purpose of this section is to describe in more detail the computation of 

S^^^fix, v) = Pe^^fix, v) = Pf(x- ^z;, v) 
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and 

Without loss of generality let us consider a translation of the form f (x — Tg{v), v). 
In addition, we fix the cell of interest as [0, h] x [0, h]. Now we are primarily interested 
in an interval of length h and thus define Pi{x) = pi{j^ — 1). Then we have 

h 

Pi{x)Pj{x)dx = ^i^Sij. 



We have to first translate and then project a function that can be expanded as 

M N 

f{x,v) = ^ ^bjnnPrn{x)Pn{v) 
m—O 71—0 

onto the finite dimensional approximation space. Our goal is to compute the coeffi- 
cients of f{x — Tg{v)^v). These are given by 

o-i3^^———Tp^——^l I Pi {x)Pj (v) f {x - Tg{v), v)dxdv 



{2U 


Kl)(2j + 1) 








Kl)(2j + 1) 




h 




Kl)(2j + 1) 


h 



JO 



where 



2_^b,nn I Pj{v)Pn{v) - I Pl{x) 
^&mn / Pj{v)Pn{v)Hira{g{v)T/h)dv, 



Him{5) ^ \ t Pii^)Pm{x - 5h) Ax, 6 = 



Tg{v)) dx dv 



(2.7) 



For a fixed v the function iJ;,„ can be evaluated explicitly. This is done, up to 
order 3, in p?7j- We will instead use a Mathematica program which can generate a 
representation of Hi^ (up to arbitrary order) in C code that can then be embedded in 
the C++ implementation. Note that it is sufficient to only evaluate Hi,n for < S < 1 
as the negative values of 6 follow by a symmetry argument and integer multiplies 
correspond to a shift of the cells only. Also, the computation of Him{5) for — 1 < ^ < 1 



shows that only two terms from the sum in (2.71 do not vanish. That is, we need only 
the data from the same cell as well as a single neighboring cell (either the right or left 
neighbor) to compute an application of a splitting operator. This follows easily from 
the fact that the support of the Legendre basis functions are within a single cell only. 
More details are given in Tf]. 



It remains to evaluate the remaining integral in equation (2.7). Since g{v) is at 
most a polynomial of degree € (in a single cell) we have to integrate a polynomial of 
degree at most f.'^. We use a Gauss-Legendre quadrature rule of appropriate order. 

Note that in order to guarantee the stability of our scheme it is of vital importance 



that we can compute the exact result of the integral in equation (2.7). If only an 



approximation is used instabilities can occur (see section 3.4 and |18) ) 
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2.7. Polynomial approximation of functions with a small jump discon- 
tinuity. In this section our goal is to prove a bound concerning the approximation 
of piecewise polynomials of degree t with a single jump discontinuity. For notational 
simplicity we will be concerned with a function of a single variable only; the general 
case is a simple tensor product of the situation described in this section. Thus, the 
operator P is here understood as the orthogonal projection with respect to the one- 
dimensional Legendre polynomials of degree less or equal to I. The starting point of 
our investigation is the result in Theorem |2.3[ which is applicable only if we can as- 
sume that (7 is ^-|- 1 times continuously differentiable. This assumption is not satisfied 
for the discontinuous Galerkin approximation considered in this paper. However, we 
will use the result as a stepping stone to prove a similar bound for the approximation 
of functions with a small jump discontinuity. 

Theorem 2.3. Suppose that g e C^+i([0, h]). Then 



for all fc e {0... ,^}. 

Proof. In ^91 p. 59] it is shown that Pg — g changes sign £ + 1 times. From 
this, it follows that {PgY'^^ — g'-'^'-' changes sign i + 1 — k times. Therefore, {Pg)^^^ 
is an interpolation polynomial of g^''^ of degree i + 1 ~ k. Using the standard error 
representation for polynomial interpolation we get the desired result.D 

For numerical methods that rely on a smooth approximation of the solution (for 
example, using Hermite or spline interpolation as in [5]) sufficient regularity in the 
initial condition implies the bound given in Theorem |2.3| for any approximation that 
has to be made in the course of the algorithm. 

This assumption, however, is violated if we consider a discontinuous Galerkin ap- 
proximation as, even if the initial condition is sufficiently smooth, the approximation 
will include a jump discontinuity at the cell boundary. Thus, we are interested in a 
bound that still gives us an equivalent result to that stated in Theorem |2.3| in the 
case of a function with a small jump discontinuity. The following theorem is thus the 
central result of this section. For simplicity we consider a single cell only. 

Theorem 2.4. Suppose that g: [Q,h] — > M is piecewise polynomial of degree 
I with a single discontinuity at xq G [0, /i]. In addition, we assume that the jump 
heights = g^''^ixo+) - g^^\xo-) satisfy leC^^I < ch^-^+^ for all k G {0,...,^}. 
Then, 



gi^) _ (pg)(^) 



< Ch'~^+\ 



for all k S {0, . . . Note that the constant C only depends on c and the constant 
in Theorem \2.9\ 

Proof. Let us assume that xo G (0, /i) (otherwise the result is immediate). We 
smooth the piecewise constant function g^^^ in the following way 

p^'\x) = '^x + g^'\Q). 



Now, upon integration we get 

h (£+1) 



fe=0 
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where we choose the coefRcients in such a way that the Taylor polynomial of g ex- 
panded at ma 
representation 



panded at matches the first £ terms of g, i.e. ak — ^"ir^- This gives us the following 



(2.8) 



Now let us consider the integral (for x > Xq) 

+ g^"'''\xo+) - g^^-'\xo~~) 
= p("-i) (x) - (x) + e("-i) , 

where the last identity follows from the choice we made above. Now we know that 

(for S£_i > Xq) 

" P^'\se) - g^'Ks,) ds, ^ p(^-i)(.,„i) - g^'^'\s,^,) + eC^'D 
and further (for Sf„2 > Xq) 

By an induction argument we can then estimate the approximation error as 



\p{x) -g{x)\ < 



^0 



p^'\sf)~g('\s,)Asi...ds2ds^ 



< 

/o Jo 







fc=0 



P*^^ (se) - .g'^' (s£) ds^ . . . dsadsi + ch^^^ 



In addition we easily follow from equation (2.8) that 



< c. 



Now let us estimate the approximation error 



^'5-511 = WPg - Pp + Pp - p + p - gW 

< \\P{g-p)\\ + \\Pp-p\\ + \\p-g\\ 

< Ch'+\ 



where in the last line we have used Theorem 12.31 and the well known fact that the 
projection operator P is a bounded operator in the infinity norm. The latter can be 
seen, for example, by estimating (2.5). 
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To get the corresponding resuh for the fcth derivative we foUow largely the same 
argument. The last estimate is then given by 



{Pg) 



(k) 



<Ch-'^\\P{g^p)\\+Ch 



(fe) _ 



9 

•-k+l 



(Pp) 



(fe) 



P 



(fe) 



where the estimate for the first term follows by the well-known Markov inequality (see 

e.g. m)- □ 

Let us discuss the principle of applying Theorem |2.4| First the operator P is 
applied to f{jT), i.e. a point on the exact solution, and we can assume the necessary 
regularity to apply Theorem |2.3[ Consequently, we get a jump discontinuity of heights 
at most 



<2 /«(jt)-P/WOt) 



< Ch 



!-fe+l 



< Ch'^-''+\ 0<k<£. 



Now the projected function is translated by a splitting operator (the example g{x) = 
{Pf{jT)){x — vt) is illustrated in Figure 2.1 1 and projected back on the finite dimen- 
sional subspace. The resulting error up to the ^-th derivative is then given by (see 
Theorem [Oil 



iPg) 



<ch'- 



fe+i 



(k) 

From this we can also follow that the new jump heights e{ are at most 
le^'^^l < 2||(P.9)(''') - g'^'^'^W < Ch^-''+\ < k < i. 





cell 1 



cell 2 



cell 1 



cell 2 



Fig. 2.1. Projected smooth function with a jump discontinuity at the cell boundary (left) and 
translation with a discontinuity inside the cell (right). Only two cells in the x direction are shown. 

Since we only have to repeat this procedure a finite number of times (i.e. for a 
single step of the Strang splitting algorithm) and the assumptions of Theorem 2.4 are 



satisfied uniformly for all /(t), we can find a uniform constant C such that the desired 
estimate holds. 

Strictly speaking this argument is only valid for a constant advection (i.e. where 
V is fixed). However, we can always decompose the projection operator as P = PyPx] 
that is, into a projection in the x-dircction (that depends on the parameter v) and 
a subsequent projection in the w-direction. Due to the special form of the advection 



we consider (see section 2.61, the projection in the x-direction gives a function that is 



piecewise polynomial in every cell. Thus, the projection in the u-direction poses no 
difficulty. 
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2.8. Consistency. It is the purpose of this section to formulate assumptions 
under which we can show a consistency bound for the initial value problem given in 
equation (2.1 1. For notational simplicity, we will denote in this section the solution 
of (2.1 1 at a fixed time tk = kr by /q. The function /o defined as Pfo is a (possible) 
initial value for a single time step (i.e., a single application of the splitting operator 
Sk or Sk)- Since we consider consistency we are interested in the non-linear operator 
S that is given by 

5(.) = 5(^)e-^(^*"'(-))5(^)(.), 
and the corresponding spatially discretized operator 

^(•)-^(-^)Pe^^(^'"'(-))^(^)(.). 

Let us first give a simple consequence of the variation-of-constants formula. 
Lemma 2.5. Suppose that fo is absolutely continuous and that /i/2,/i/2 
integrable. Then 



Proof. Let ^(r) = S^^^^^/^^^h. Then 

g' - B{h/^)g = BCh/2)9 + (s(/i/2) - BCh/2)) 5, 
which can be rewritten by the variation-of-constants formulas as 



from which the desired result follows immediately. □ 

The next two lemmas will be the crucial step to prove consistency. First, we con- 
sider the error made by the (exact) splitting operators due to the space discretization. 
Note that the assumptions are exactly the same as those needed for the development 
in section [221 to hold. 

Lemma 2.6. Suppose that Jq G Cp+.^j,. Then 

\\Sfo-Sfo\\<C{Th'+' + h'+'), 

where C depends on \\fo\\Qe+i- (but not on r and h). 

Proof. Let us define /1/2 = S^^^/q. Then, we can write 



< S'('^) (^5(s(/i/2)) _ 5'(B(A/2))^ 5'(^) 
By using Lemma |2.5| and the definition of B we get 

< CT\\E{fy,) - i?(A/2)|| max (e'^^^f^^-^ S'^^^ fo 
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Note that the set of measure zero, where dy is not defined, does not influence the 
estimate in the L°° norm as we are only concerned with equivalence classes of (es- 
sentially) bounded functions. Finally, since E is given by equation (2.3) it follows 
that 

\\E{fi/2) - E{h;,)\\ < 11/1/2 - A/2II < ll/o - /oil < Ch'+\ 

which concludes the proof. □ 

Second, we consider the error made due to the approximation of the (exact) 
splitting operators. Note that the assumptions are exactly the same as those needed 



for the development in section 2.7 



Lemma 2.7. Suppose that fo e Cp+.^(,. Then 

\\Sfo-Sfo\\<C{Th'+' + h'+'), 

where C depends on ||/o||^*+i (hut not on r and h). 
Proof. Let /1/2 — S'^^'f^. Then, we can write 



Now we estimate the three terms independently. The estimation of the first and third 
term is straightforward. We get 



(P-l) (5(^)5(^(/v2))5(^)/g 



and 



S(A)siB(h,2)) (^S^A) _ ^(A)^j 



<c||(5(^)-^(-^))/o 

= cj{l~P)S^^Uo 
< Ch'+\ 

To estimate the second term we employ Lemma |2.5| which gives 

SiA) ^^^(^(/i/s)) „ ^(S(/i/2))^ 5(A) 



< CT\\E{fy,) ~ E{f,/2)\\ max 



do, 



Note that the set of measure zero, where is not defined, does not influence the 
estimate in the L°° norm as we are only concerned with equivalence classes of (essen- 
tially) bounded functions. As in the last lemma E is given by equation (2.3) and thus 
it follows that 



\E{f- 



1/2) 



E{f: 



1/2) 



<ll/: 



1/2 



/: 



1/2I 



il-P)fi/2\\<Ch 
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which concludes the proof. □ 

Theorem 2.8 (Consistency). Suppose that fo e Cpc^^^^~^^'^\ Then 

\\Pf{h)-Sfo\\<C{r' + Th'+'+h'^+'), 

where C depends on ||/o|||^max(*+i,3) (but not on t and h). 
Proof. We write 

\\Pf{h) ~ ~Sh\\ = \\Pf{h) - PSfo + PSfo - 5/0 + Sfo - Sfo + Sfo - Sf4 

< \\P{f{h) - S.M II + ||P5/o - 5/0II + \\Sfo - 5/0II + ||5/o - SfoW 

< Ct^ + Ch'+^ + ||5/o - SM + l|5/o - ~SM, 

where the first term was bounded by Theorem 4.9 in The two remaining terms 
can be bounded by Lemmas 2.6 and 2.7 to give the desired estimate. □ 

2.9. Convergence. To show consistency it was most convenient to bound all 
terms in the infinity norm. Bounds in the or norms then follow since we 
consider a compact domain in space and velocity. However, for stability (and thus 
convergence) we need to bound the operator norm of the projection operator P by 1. 
Since such a bound is readily available in the norm (as an orthogonal projection is 
always non-expansive in the corresponding norm) we will use it to show convergence. 
Note that this is not a peculiarity of our discontinuous Galerkin scheme. For example, 
in [S] stability for two schemes based respectively on spline and Hermite interpolation 
is shown in the norm only. 

Theorem 2.9 (Convergence). For the numerical solution of (2.1 1 we employ 
the scheme (2.6). Suppose that the initial value fo G (jmax{£+i,3} rion-negative and 
compactly supported in velocity. Then, the global error satisfies the bound 



sup 

0<n<JV 



finr) 



\k=0 



< C 



+ h'+' 



where C depends on T hut is independent of t, h,n for < nr < Nt = T. 
Proof. From (2.6) we get 



Now we can derive a recursion for the error in the norm 

e«+i = ll/n+i - f{nT + t)||2 

< ||5(^)pe-^(/" + V2)^(^) _ S^A)p^rB(P.i-pfinr))giA) j^^^y 

+ ||^(-4)pc-^(^"*-'^/("^))^(^)(/„ - /(nT))||2 

+ ||^(-4)pe-^(^"*-'^/("^))^(^)(l - P)f{nT)h 

+ ||^(-4)Fc-^(^"*''^/("^))5(^)p/(nT) -Pf{nT + r)||2 
+ ||(P-l)/(nr + T)||2. 
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The first term is estimated with the help of Lemma |2.5[ the fourth one with Theo- 
rem [2]8) This gives us 



en+i < CT\\fn+i/2 - Pe 5 ^F/(nr) 1 1 2 + ||/„ - /(nr)|l2 + C (r^ + Th'+' + h'+^] 
< (1 + CT)ek + C{t' + Th'+' + h'+^) . 

Applying a discrete Gronwall lemma to the above recursion then gives 

<C t' + h'+' + — 



which is the desired bound as the constant C can be chosen uniformly in [0,T]. This 



follows from the regularity result (Theorem 2.1 1 which gives us the desired bound 



for Theorem 2.8 if fo € Qma.x{e+i,3} non-negative and compactly supported with 
respect to velocity. □ 

2.10. Extension to higher dimensions. In three dimensions the splitting 
scheme is given by (for simplicity we consider a single time step only and thus drop 
the corresponding indices) 

S^^^f{x,v)^f(x-^v,vy (2.9) 
^(^)/(a;, v) - f{v, X - TE{f,/2, x)). (2.10) 



The expression in equation ( |2.9[ ) can be easily decomposed into three translation in a 
single dimension, i.e. 

with Ax = -Vxdx, Ay = -Vydy, and A^ = -Vzdz- 



The discussion is more subtle for the expression in equation (2.101. In this case 
we can still use the decomposition given above; however, if we introduce a space 
discretization we will have to project back not onto a 1 -|- 1 dimensional space but 
onto 1-1-3 dimensional space. This is an important implementation detail; however, 
the convergence proof is (except for notational difficulties) unaffected. 

As most of the derivation in this paper is conducted within the framework of 
abstract operators, the extension to multiple dimensions is straightforward. In Lem- 
mas 



2.6 and |2.7| we have to consider a more general differentiation operator (i.e. a 
directional derivative). However, this represents no difficulty as the existence and 
regularity results are not restricted to the 1-1-1 dimensional case (see [T^). 

Therefore, it remains to generalize the discussion given in [9] to multiple dimen- 
sions in the case of the Vlasov-Poisson equation. The abstract results hold inde- 
pendently of the dimension and the specific details of the operators A and B. The 
remaining computations are somewhat tedious, however, as the existence and regu- 
larity results are essentially the same the proof can be extended in a straightforward 
fashion. 

3. Numerical simulations. The purpose of this section is to perform a number 
of numerical simulations in order to establish the validity of the implementation. 
The recurrence phenomenon in the context of higher order implementations in space 
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is discussed in section |3.1| In section |3.2| the order of the method in the strong 
Landau damping problem is investigated. We will also reproduce some medium time 



integration results for linear Landau damping (section 3.3 ) and investigate the stability 



for the Molenkamp-Crowley test problem (section 3.4 1. 

The computer program used for the numerical simulations performed in this sec- 
tion is implemented in C++. It employs heavily the concept of templates and operator 
overloading to provide a compact implementation that is easily extendable to the 
multi dimensional case. As a result, the core program consists of only about 800 lines 
of source code (excluding unit tests but including all the logic needed to carry out 
the simulations in this section) while still maintaining an implementation with good 
performance. 

3.1. Recurrence. It is well known that piecewise constant approximations in 
velocity space lead to a recurrence phenomenon that is purely numerical in origin. This 
behavior has been investigated, for example, in and [T5] . In [53] it is demonstrated 
by a number of numerical experiments that in the weak Landau damping problem 
this phenomenon is also purely a numerical artefact. 

From an analytical point of view the recurrence phenomenon is most easily un- 
derstood for an advection equation, i.e. a function /(t, x, v) satisfying the equation 

dtf - -vd^f. (3.1) 

For its numerical solution consider a piecewise constant approximation of / in ve- 
locity space. This approximation results in slices in velocity space that correspond 
to the average velocity in a particular cell. Let us further assume that the velocity 
space [— Wmax, Wmax] cousists of an odd number of cells and that the interval [0,47r] is 



employed in the space direction. Then the solution of (3.1) is a periodic function in 
time and the period p is easily determined to be 

hyp = An. 

That this is only a numerical artefact is a simple consequence of the fact that p tends 
to infinity as hy tends to 0. However, for the purpose of this section it is instructive 
to compute the exact solution for the following initial value 

foix, v) = (l + 0.01 cos(0.5a;)) . 



The solution of (3.1) is then given by 



e-" /2 

f{t, X, v) = (1 + 0.01 cos(0.5a; - 0.5vt)) . 

V Stt 

This function, however, is not periodic in time (with a period being independent of 
v). To represent this more clearly, we compute the electric energy 

Sit) = E{t,xrdx = ^e-0■2^*^ (3.2) 

where the electric field E{t,x) is determined as before by 

E{t, x)= I K{x, y)( I f{t, y, v)dv ~ l] dy. 
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Note that the kernel K{x,y) is defined in equation (2.3). Thus, the electric energy is 
exponentially decreasing for the exact solution (but periodic in time for the numerical 
solution). One might naively expect that this phenomenon vanishes as soon as one 
considers an approximation of degree at least 1 in the velocity direction. While it is 
true that the solution is no longer periodic, as can be seen from Figure 3.1 errors 



in the velocity approximation still result in a damped recurrence of the electric field. 
Note that the size of this recurrence effect seems to be determined by the space 
discretization error. 



o 
c 



<0 

u 



o 
o 

i— 
\— 
01 
(D 
> 

to 

0) 



0.001 



0.0001 



Nv=32, 1=0 - - ■ 

Nv=64, 1=0 

Nv=32, 1=2 

exact solution 




Fig. 3.1. Recurrence phenomenon for the advection equation (top). Note that while there is no 
periodicity in the second order approximation a recurrence-like effect from the finite cell size is still 
visible. The (absolute) error as compared with the exact solution given in | |3.2| in the discrete 
norm (bottom) behaves as expected. In all simulations 32 cells and an approximation of order 2 (i.e. 
£ = 1) have been employed in the space direction. The number of cells and the order of discretization 
in the velocity direction is indicated in the legend. In all computations t = 0.05 is used. 
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Fig. 3.2. The recurrence phenomenon for Landau damping. Note that even though a higher 
order of approximation in the velocity direction improves the solution a recurrence- like effect is still 
visible. In all computations r = 0.05 is used. 



3.2. Order. We can investigate both, the order of convergence in time (i.e. where 
the space error is small enough over the range of step sizes r we are interested in) and 
the order of convergence in space (i.e. where the step size is chosen small enough such 
that the time integration error is negligible). The order of the time integration has 
already been investigated in [9]. Thus, we focus on the convergence order in space. 

Let us consider the Vlasov-Poisson equations in 1+1 dimensions together with 
the initial value 

fo{x,v) ^ — Le"^'/2(l + acos(0.5x)). 

This problem is called Landau damping. For a = 0.01 the problem is called linear 
or weak Landau damping and for a = 0.5 it is referred to as strong or non-linear 
Landau damping. As can be seen, for example, in [51lllj and |21j Landau damping is 
a popular test problem for Vlasov codes. 

In our numerical simulations, all errors are computed with respect to a reference 
solution (such as to exclude unwanted effects from the time discretization). The 
reference solution uses 512 cells with £ = 2 and a time step of t = 0.1. The results 



for strong Landau damping are given, up to order 3, in Figure 3.3 It can be seen 
that the accuracy improves with the desired order as the cell size decreases. Thus, 
the results are in good agreement with the theory developed in this paper. 

3.3. Landau damping. The Landau damping problem has already been intro- 
duced in the previous section. In this section we are not interested in the desired 
order of the numerical algorithm but in the comparison with the exact solution of the 
Vlasov-Poisson equation. However, since an exact solution of the full Vlasov-Poisson 
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Fig. 3.3. Error of the particle density function /(I,-,-) for non-linear Landau damping on 
the domain [0, Att] x [—6, 6] as a function of N , the number of cells in both the space and velocity 
direction. 



equation is not known we will instead use a result that gives us the asymptotic decay 
rate 7 of the electric field in the case of weak Landau damping (see e.g. [T]). Thus, 
we compare the decay of the energy stored in the electric field with the graph of 
e~27t^ where 7 sa 0.1533. A number of such simulations have already been conducted 
(see e.g. [13]). However, due to the recurrence eff'ect usually a large number of cells 
have to be employed in order to get accurate results for medium to long time inter- 
vals. For reference we note that [53] uses — — 1024 whereas [13] uses up to 
Nx — 2000, A'^ — 1600. The results of our simulation are shown in Figure [3^ (the 
number of cells is ^ ^ 256 for £ = 1 and A^^ = A^,. = 128 for £ ^ 2). This 
experiment clearly shows that high-order approximations in space and velocity pay 
off. 

3.4. Stability. In advection dominated problems instabilities can occur if the 
numerical integration is not performed exactly. In [IB] this is shown for the Molenkamp- 
Crowley test problem, i.e. 



dtfit, X, y) = 2TT{ydx - xdy)J{t, x, y) 
fiO,x,v) = faix,y), 



(3.3) 



where 



fo{x,y) 



cos^(27rr) t < \ 
else 



with = [x -\- \Y + y'^ . The solution travels along a circle with period 1. We will 
solve the same problem using the algorithm presented in this paper and show that no 
instabilities occur if a quadrature rule of appropriate order is used. This results are 



given in Figure 3.5 Note that this is exactly what is expected based on the theoretical 
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Fig. 3.4. The decay of the electric field is shown for = Ny = 256, £ = 1 (top) and 
Nx = Ny = 128, 1 = 2 (bottom). In both cases a relatively large time step of r = 0.2 is employed. 



analysis done in section [2] However, it is not true that such stabiHty results hold for 
arbitrary schemes (see e.g. |!18! , where a finite element scheme of order 2 is shown to 
be unstable for most quadrature rules). 

4. Conclusion. In the present paper we have extended the convergence analysis 
conducted in 9J to the fully discretized case using a discontinuous Galerkin approx- 
imation in space. The results are only presented in case of the 1 + 1 dimensional 
Vlasov-Poisson equation. However, we have given a short argument that the ex- 
tension to multiple dimensions is although tedious in principle straightforward. In 
addition, we have presented a number of numerical simulations that investigate the 
behavior of the proposed algorithm. These simulations suggest that the algorithm 
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Fig. 3.5. Stability for the Molenkamp Crvwlay test problem. The initial value is displayed at 
the top. The numerical solution after 60 revolutions with t = 0.02, A^^ = Ny = 40, and £ = 2 is 
shown at the bottom. As expected no numerical instabilities are observed. The negative values in the 
numerical solution are a consequence of the space discretization error and are propagated in space 
by the numerical algorithm (see the complex contour line of at the bottom). However, this fact 
has no influence on the stability of the scheme. 
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has certain advantages over similar algorithms that employ a piecewise constant ap- 
proximation in space and velocity. 
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